#!/usr/bin/perl
# 
# Exercise 8.12
# 
# Write a program that, given two amino acids, returns the probability that a single 
# mutation in their underlying (but unspecified) codons results in the codon of one 
# amino acid mutating to the codon of the other amino acid.
#
# Answer for Exercise 8.12

# Plan:
# Get all the codons that code for the first amino acid.  Then, for each codon,
# there are nine possible mutations.  Simply check each one to see if it results
# in a codon for the second amino acid.  Finally, divide the number of codons
# that you find do code for the second amino acid, by the total number of mutated codons
# considered (9, 18, 27, or whatever), to get the probability.

use strict;
use warnings;
use BeginPerlBioinfo;

srand(time|$$);

print "Give a one-character amino acid: ";
my $aa1 = <STDIN>;
chomp $aa1;

print "Give a different one-character amino acid: ";
my $aa2 = <STDIN>;
chomp $aa2;

my $probability = mutation_probability($aa1, $aa2);

print "The probability that amino acid $aa1 mutates to amino acid $aa2 is $probability\n";

exit;

sub mutation_probability {

	my($aa1, $aa2) = @_;

	# To count the mutations of any codon for the first amino acid, that results in a codon for
	# the second amino acid.
	my $count = 0;

	my(%reverse_genetic_code) = (
	
	'A' => 'GCA GCC GCG GCT',		# Alanine
	'C' => 'TGC TGT',			# Cysteine
	'D' => 'GAC GAT',			# Aspartic Acid
	'E' => 'GAA GAG',			# Glutamic Acid
	'F' => 'TTC TTT',			# Phenylalanine
	'G' => 'GGA GGC GGG GGT',		# Glycine
	'H' => 'CAC CAT',			# Histidine
	'I' => 'ATA ATC ATT',			# Isoleucine
	'K' => 'AAA AAG',			# Lysine
	'L' => 'CTA CTC CTG CTT TTA TTG',	# Leucine
	'M' => 'ATG',				# Methionine
	'N' => 'AAC AAT',			# Asparagine
	'P' => 'CCA CCC CCG CCT',		# Proline
	'Q' => 'CAA CAG',			# Glutamine
	'R' => 'CGA CGC CGG CGT AGA AGG',	# Arginine
	'S' =>'TCA TCC TCG TCT AGC AGT',	# Serine
	'T' => 'ACA ACC ACG ACT',		# Threonine
	'V' => 'GTA GTC GTG GTT',		# Valine
	'W' => 'TGG',				# Tryptophan
	'Y' => 'TAC TAT',			# Tyrosine
	'_' => 'TAA TAG TGA',			# Stop
	);

	# The set of codons that code the first amino acid
	my @codons = split(' ', $reverse_genetic_code{$aa1});

	foreach my $codon (@codons) {

		# There are nine such mutations for each codon
		foreach my $mutation (mutate_codons($codon)) {

			if (codon2aa($mutation) eq $aa2) {

				$count++;
			}
		}
	}

	# return the probability
	return $count / (9 * @codons);
}


#
# A helper subroutine: given a codon, it returns all the codons that
# can result from a single base change -- 9 in all.
#

sub mutate_codons {

	my($codon) = @_;

	my @bases = ('A','C','G','T');

	# To collect the output mutations
	my @mutations = ();
	
	#
	# Note the loop within a loop.  The outer loop loops through the 3 positions in the
	# codon.  The inner loop loops over the possible bases.
	#
	for(my $i = 0; $i < 3 ; ++$i) {
		foreach my $base (@bases) {

			my $mutation = $codon;

			if ($base eq substr($codon,$i,1)) {
				next;
			}else{
				substr($mutation,$i,1) = $base;
			}
			push(@mutations, $mutation)
		}
	}
	return @mutations;
}
